subroutine rhs(u,beta,u_rhs)
implicit none
real*8::u(:,:,:),u_rhs(:,:,:),beta(:,:,:,:),one=1.
integer::ix,iy,iz,nx,ny,nz,sx,sy,sz
real*8,parameter::c=1.,dx=1.,dy=1,dz=1

nx=size(u,1)
ny=size(u,2)
nz=size(u,3)

do ix= 4, nx-3
  do iy= 4, ny-3
    do iz= 4, nz-3
              sx=sign(one,beta(ix,iy,iz,1))
              sy=sign(one,beta(ix,iy,iz,2))
              sz=sign(one,beta(ix,iy,iz,3))
              u_rhs(ix,iy,iz)=-sx*1./12.*beta(ix,iy,iz,1)*(3.*u(ix+sx,iy,iz)+10.*u(ix,iy,iz)&
                                -18.*u(ix-sx,iy,iz)+6.*u(ix-2*sx,iy,iz)-u(ix-3*sx,iy,iz))/dx&
                             -sy*1./12.*beta(ix,iy,iz,2)*(3.*u(ix,iy+sy,iz)+10.*u(ix,iy,iz)&
                                -18.*u(ix,iy-sy,iz)+6.*u(ix,iy-2*sy,iz)-u(ix,iy-3*sy,iz))/dy&
                             -sz*1./12.*beta(ix,iy,iz,3)*(3.*u(ix,iy,iz+sz)+10.*u(ix,iy,iz)&
                                -18.*u(ix,iy,iz-sz)+6.*u(ix,iy,iz-2*sz)-u(ix,iy,iz-3*sz))/dz
    end do
  end do
end do

end subroutine



subroutine RungeKutta(u,n,beta,uNext)
implicit none
integer,intent(in)::n
real*8::u(:,:,:),beta(:,:,:,:),uNext(:,:,:),u_k1(n,n,n),u_k2(n,n,n),u_k3(n,n,n),u_k4(n,n,n),u_tmp(n,n,n)
real,parameter::dt=0.25
interface

subroutine rhs(u,beta,u_rhs)
implicit none
real*8::u(:,:,:),u_rhs(:,:,:),beta(:,:,:,:)
integer::ix,iy,iz,nx,ny,nz
real*8,parameter::c=1.,dx=1.,dy=1,dz=1
end subroutine

end interface


 call rhs(u,beta,u_k1)
 
    u_tmp=u+0.5*dt*u_k1 

  call rhs(u_tmp,beta,u_k2)
    u_tmp=u+0.5*dt*u_k2

  call rhs(u_tmp,beta,u_k3)
    u_tmp=u+dt*u_k3

  call rhs(u_tmp,beta,u_k4)

     uNext=u+(1./6.)*dt*(u_k1+2*u_k2+2*u_k3+u_k4)

  end subroutine


program main
implicit none
integer,parameter::nx=20,ny=20,nz=20,m=50
real*8:: u(m,nx,ny,nz)=0,u_k1(nx,ny,nz)=0,u_k2(nx,ny,nz)=0,u_k3(nx,ny,nz)=0,u_k4(nx,ny,nz)=0,u_tmp(nx,ny,nz)=0,beta(nx,ny,nz,3)=0
integer::it,ix,iy,iz

interface

subroutine RungeKutta(u,n,beta,uNext)
implicit none
integer,intent(in)::n
real*8::u(:,:,:),beta(:,:,:,:),uNext(:,:,:),u_k1(n,n,n),u_k2(n,n,n),u_k3(n,n,n),u_k4(n,n,n),u_tmp(n,n,n)
real,parameter::dt=0.25
end subroutine


end interface
beta(:,:,:,1)=1
beta(:,:,:,2)=1

do iy=16,20
do ix=1,iy-1
  u(1,ix,iy-ix,:)=1
end do
end do

do it=1,m-1
  call RungeKutta(u(it,:,:,:),nx,beta,u(it+1,:,:,:))
end do

open(file='a.txt',unit=10,action="write")
do it=1,m

 write(10,"(20f9.5)")u(it,:,:,10)
 write(10, * )

end do

open(file='b.txt',unit=11,action="write")
do it=1,m

 write(11,"(20f9.5)")u(it,:,10,10)

end do
end program
